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Linsley used the Volcano Ranch array to collect data on the lateral distribution of showers produced 
by cosmic rays at energies above 10^^ eV. Very precise measurements of the steepness of the lateral 
distribution function were made on 366 events. The current availability of sophisticated hadronic 
interaction models has prompted an interpretation of the measurements. In this analysis we use 
the AIRES Monte Carlo code to generate showers, together with GEANt4 to simulate the detector 
response to ground particles. The results show that, with the assumption of a bi-modal proton and 
iron mix, iron is the dominant component of cosmic rays between 5 x 10^^ and 10^^ eV, assuming 
that hadronic interactions are well-described by QGSJET at this energy range. 
PACS: 96.40, 13.85.T High Energy Cosmic rays; mass composition 

I. INTRODUCTION 

The measurement of the mass composition of cosmic rays above 10"'^^ eV is a challenging problem. This information 
is as important as the energy spectrum and the anisotropy in determining cosmic ray origin. One must know the likely 
mass range of a particular data set before one can interpret anisotropy information confidently, given the influence 
of galactic and intergalactic magnetic fields. Our knowledge of the mass composition of cosmic rays above 10^^ eV 
remains very limited. Recent re-interpretation of measurements of the lateral distribution of water-Cerenkov signals 
made at Haverah Park [1] suggests a composition of 34% protons and 66% iron in the range 2 x 10"'^^ - 10^^ eV. This 
contrasts with earlier claims, from observations made using Fly's Eye, that the composition changes from a heavy 
mix around 3 x 10^^ eV to a proton dominated flux around 10^^ eV [2]. At Yakutsk, both the inferred values of the 
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depth of shower maximum {X^ax) and the muon density favor a composition change from a mixture of heavy and 
Ught components to Ught composition over the same energy region [9]. Prom HiRes/MIA data [7], there are claims 
that there is a rapid change from a heavy to a hght composition between 0.1 and 1.0 EeV. A recent analysis by the 
HiRes collaboration of data collected in the energy range between 10^^ - 10^^ '* eV [8] is consistent with a nearly 
constant purely protonic composition. The fraction of protons, however, decreases when they interpret their data 
using SIBYLL2.1. On the other hand, the AGASA group have argued for a "mixed" unchanging composition from 
1-10 EeV [4] (using MOCCA for the simulations). A recent analysis of the muon component in air showers with 
aires/qgsjet98 around 10^^ eV by the AGASA collaboration indicates a relatively light average composition [5]. 

The source of the discrepancy between different experiments is not understood and it is important to resolve the 
issue, because of its implications for cosmic ray models of origin, acceleration and propagation. Volcano Ranch data 
may provide a path for further understanding. 

Following the successful re-examination of the Haverah Park data [1] with modern shower models, we report here 
a similar analysis using the Volcano Ranch data, collected by Linsley [16] to determine the shape of the lateral 
distribution of air showers. This is the first attempt to examine the Volcano Ranch data with the results of Monte 
Carlo calculations, using Monte Carlo tools that were unavailable when the data were recorded in 1970. It is timely 
as the situation on mass composition above 10^'' eV remains confused and the steepness of the lateral distribution is 
sensitive to the depth of maximum of the shower, and therefore to the primary composition and to the character of 
the initial hadronic interactions. 

To simulate the development of the air showers, we have used the aires [10] code (version 2.4.0), with the hadronic 
interaction generator QGSJET98 [11]. The results of the simulated showers were convolved with a simulation of the 
detector response made using GEANt4 [15]. A comparison of two hadronic generators (qgsjet98 and SIBYLl2.1 ) 
was presented in [26]. Both give satisfactory descriptions of the data, but we have preferred to use qgsjet98 because 
this model has been shown to be consistent with experimental data at energies up to 10 PeV and beyond [13,14] 

II. THE VOLCANO RANCH ARRAY 

The pioneering Volcano Ranch instrument consisted of an array of scintillation counters. The array was operated 
in three configurations from 1959-1976 at the MIT Volcano Ranch station located near Albuquerque, New Mexico 
(atmospheric depth 834 g cm~^). One of its many distinctions was the detection of the first cosmic ray with an 
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energy estimated at lO^^eV [19]. The final configuration, of relevance here, comprised 80 detectors of surface area 
0.815 m^, scintillator thickness of 9.032 g cm~^ laid out on a hexagonal grid with a separation of 147 m (Figl). This 
configuration allowed precise measurement of the lateral distribution of the detector signals. The steepness of the 
lateral distribution, and its fluctuations, can be used to explore the primary mass composition as in [1]. Fortunately, 
in his various writings, Linsley has left unusually detailed descriptions of his equipment, together with examples of 
events and a description of his data reduction methods. 
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FIG. 1. Volcano Ranch array in the final configuration. 



III. LATERAL DISTRIBUTION FUNCTION 

A generalized version of the Nishimura-Kamata-Greisen (NKG) formula was used to describe the lateral distribution 
of particles at ground in minimum ionizing particles per square metre (mipsm^^) for Volcano Ranch data [24]. This 
lateral distribution function is given as 



Nfit s / r \ I ^ T 



-(»?-") 



5vR(r) = ^C(a,r/) — 1 + — (1) 

normalized to shower size A^/it with 

^ r(r? - a) 

27rr(2 - a)r(77 - 2) ' ^ ' 

Here Vm is the Moliere radius, which is 100 m for the Volcano Ranch elevation. r\ and a are parameters that 
describe the logarithmic slope of this function. 



From a subset of 366 showers detected with the array, the form of 77 as a function of zenith angle 6 and shower size 
Nfit was found to be [20] : 

Nfit)) =a + b{sece - 1) + clogio{^) (3) 
with a = 3.88 ± 0.054, b = -0.64 ± 0.07, and c = 0.07 ± 0.03 where a fixed value of a = 1 was adopted. 

IV. SIMULATION OF THE DETECTOR RESPONSE OF THE VOLCANO RANCH ARRAY 

The AIRES code provides a realistic air shower simulation system, which includes electromagnetic algorithms [18] 
and links to different hadronic interactions models. As mentioned above, we have used the QGSJET98 model for 
nuclear fragmentation and inelastic collisions. For the highest energy showers, the number of secondaries becomes so 
large that it is prohibitive in computing time and disk space to follow and store all of them. Hillas [17] introduced 
a non-uniform statistical sampling mechanism which allows reconstruction of the whole extensive air shower from a 
small but representative fraction of secondaries that are fully tracked. Statistical weights are assigned to the sampled 
particles to account for the energy of the discarded particles. This technique is known as "statistical thinning" . The 
AIRES code includes an extended thinning algorithm, which has been explained in detail [10]. The present work has 
been carried out using, in most cases, an effective thinning level eth = Eth/Eprim = 10~^ which is sufficient to avoid 
the generation of spurious fluctuations and to provide a statistically reliable sample of particles far from the shower 
core. All shower particles with energies above the following thresholds were tracked: 90 keV for photons, 90 keV for 
electrons and positrons, 10 MeV for muons, 60 MeV for mesons and 120 MeV for nucleons and nuclei. 

We have generated a total of 1735 proton and iron showers with zenith angles in the range secO = 1.0 - 1.5 and 
primary energies between 10^'' eV and 10^® eV, to match the Volcano Ranch data. To simulate the response of 
the detectors of the array to the ground particles, we utilized the general-purpose simulation toolkit geant4. Our 
procedure follows the prescription in [23], where the detector response to electrons, gamma, and muons is simulated 
in the energy range 0.1 to 10^ MeV and for five bins per decade of energy. The results of air shower simulations are 
convolved with the detector response to obtain the scintillator yield expressed in mips m~^. The computed lateral 
distributions of particles and the corresponding signal from the scintillators for photons, electrons and muons are 
displayed in Figure 2. 
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FIG. 2. Simulated lateral distributions of the three main shower components at ground level and its convolution with the 

detector response. 



A comparison between the lateral distribution measurements [24] and proton/iron showers simulated with 
AIRES/qgsjet98 including the scintillator response was presented previously [26]. Each simulated shower was thrown 
a maximum of 100 times on to the simulated Volcano Ranch array with random core positions in the range - 500 m 
from the array center. 

With the thinning method used in Monte Carlo shower generators, when particles reach the thinning energy just 
one of them is followed and multiplied by the corresponding weight at the end. Thus, to simulate the response of the 
detectors correctly, it is necessary to perform a smoothing of the densities of the ground particles around the position 
of each detector. All particles in a "sampling zone" around a given detector are selected and the statistical weight, 
as obtained from aires, is multiplied by the "sampling ratio" Adetector / ^sampling where AsampUng is the area of the 
"sampling zone" and Adetector is the corresponding detector area. This is equivalent to sampling particles on a larger 
area to get a realistic density around the detector position. As the densities depend mainly on the distance to the 
shower axis, the sampling area over which simulated particles are gathered is such that this ratio varies from about 
0.1 at 100 m to about 0.001 at 1 km. As a first check on the validity of our approach, the data of a single large 
event were compared with calculations [26]. Further checks between data and Monte Carlo were performed as the 
one in Figure 3. In this plot we present a comparison between lateral distribution measurements [24] and a 10^^'^ 
eV proton (left) and iron (right) shower simulated with AIRES/qgsjet98, including the scintillator response of the 
detectors in the Volcano Ranch array configuration. It has been shown [25] that the fluctuation of the density of 
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shower particles far from the core is quite small and that the density at 600 m, S(600m), depends only on primary 
energy. We normalize the showers to the value of S(600m) in order to decouple the normalization factor from the 
parameters related to the shape of the lateral distribution which change with primary mass but only slightly with 
shower size. The agreement between data and Monte Carlo is good and gives confidence in the procedures used. 
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FIG. 3. Comparison between lateral distribution measurements in a single event [24] and the simulated scintillator response 
in the configuration of VR array for 10^^'^ eV proton and iron showers. 



V. DERIVATION OF THE PRIMARY MASS COMPOSITION 

The nature of the primaries that initiate air showers is difhcult to establish from the average properties of the 
data. For example, an average property can be explained with a mass composition of a single species (A) or by an 
appropriate mixture of species. However, with the Volcano Ranch array, accurate measurements of j] were made on a 
shower-by-shower basis for fixed bins of zenith angle separated by 80 g cm~^ [21]. Thus the fluctuations of 77 can be 
used to break this degeneracy. Linsley determined the precision of each measurement of 77 and reported the average 
value of this quantity for each zenith angle bin. 

The average error in 77 from the fit made to the simulated lateral distributions {agim = 0.029), is smaller than the 
one reported by Linsley (cr = 0.072). To include within the simulation the effect of data reconstruction, we smeared 
each value of rj calculated by Monte Carlo using a Gaussian with a width chosen so that Linsley's overall uncertainties 
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in r] are reproduced i.e. (Jsmear = 0.0662 is found by quadrative subtraction of the average values of <Jsim and a. 

This is a minor correction since the measurement accuracy is so much smaUer than the intrinsic shower-to-shower 
spread (r.m.s.= 0.19). Thus, for each value of -q found from the Monte Carlo calculation, we have a corresponding 
and realistic estimate of its "experimental" uncertainty. We are thus able to make comparisons of Volcano Ranch 
data with our calculations. 

As a further check, we have calculated the variation of r/ with shower size and zenith angle with Monte Carlo and 
made comparisons with the Volcano Ranch data. The number of particles at ground level (Nfn) is obtained from a 
fit to the lateral distribution function (with a = 1) for fixed bins of zenith angle. The variation of rj with Nfu from 
the calculation has been compared with the average functional form of rj given by Equation 3. The results of this 
comparison for sec 9 — 1.0 - 1.1 and 1.3 - 1.4 can be seen for proton and iron showers in Figure 4. The error bars 
indicate the r.m.s. spread of data which is very much greater than the r.m.s. spread of the mean. The shaded band 
represents the fit to Volcano Ranch data, including the errors given in Equation(3) for a and b. 
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FIG. 4. Comparison of 77 as a function of shower size for sec9 — 1.0 — 1.1 (left) and secO 
aires/qgsjet98 
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1.3 — 1.4 (right) using 



The variation of 77 with zenith angle is shown in Figure 5 for events with shower size in the range log Nfit = 7.6 8.6 
(left) and log Nfit — 8.6 — 9.6 (right). One can see that the average form of 77 over a realistic range of mass composition, 
from proton to iron, is well represented by the simulations. The error bars represent the r.m.s. spread as before. 
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FIG. 5. Comparison of r; as a function of sec 6* for the first bin (left) and second bin (right) in logNfi 



A. Fitting VR data mass composition using finite Monte Carlo samples of different primaries. 

We can estimate the prinrary mass that describes the Volcano Ranch data, assuming a bi-modal composition, 
using a maximum likelihood fit for the best linear combination of pure iron and pure proton samples to match the 
data sample. The available data are in bins of i] [23]: the number of data points in several bins is small, so a 
minimization is inappropriate. A maximum likelihood technique assuming Poisson statistics was adopted. The 
probability of observing a particular mrmber of events di in a particular bin is given by exp"-^* ff^ /di\ where fi is the 
predicted value for the number of events in this particular bin. If we assume a bi-modal composition of proton and 
iron with fractions Ppe and Pp then fi =^ C{PFe + Pp) where C is the overall normalization factor between numbers of 
data and Monte Carlo events. Estimates of the fractions Pj are found by maximizing ln(L) = T.di ln(/i) — ln(c?.i!) — 
Our Monte Carlo samples are at least ten times larger than the data sample to avoid effects of finite Monte Carlo 
data size. 
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FIG. 6. The measured distributions of r] (data points) with histograms from Monte Carlo calculations of pure iron (left) and 
pure proton (right) with 1.0 < sec6' < 1.1, using QGSJET98. 
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FIG. 7. Composition fit from rj distribution for the first bin in 1.0 < sec^ < 1.1, using QGSJET98. The points are the VR 
data and the solid line corresponds to the result of the fit. 



In Figure 6 we compare the Monte Carlo results with the Volcano Ranch data points for near vertical showers. As 
can be seen the tail at large rj in the comparison with iron indicates that a lighter component is required to fit the 
experimental data. The best fit gives a mixture with (89 ± 5) % of iron, with a corresponding percentage of protons, 
and this distribution of r] is shown in Figure 7. One detail that Linsley did not describe is the distribution of shower 
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sizes that comprise the data set. What is recorded is that the median energy was lO""^^ eV and that the shower 
sizes are between 4 x 10^ and 6 x 10^. This corresponds to an energy range approximately between 10^^ eV and 
10^^ eV. To reproduce the data set, a differential energy spectrum of slope -3.5 was chosen for the primary energy 
spectrum. Simulations show that the whole area of the array was active above 10^^ eV. Using this spectrum and 
the condition that the median energy be 10^^ eV, the calculated threshold energy is approximately 10^^'^ eV. Below 
10^* eV the energy distribution of the data set was smoothed (always under the condition of median energy to be 
10^* eV) assuming that the effective area of the detector increases with energy from the threshold up to 10^^ eV. The 
conclusions about mass composition that we presently draw, are constrained by uncertainties in the details of the 
energy distribution of events recorded at VR. Additionally they are constrained by the hadronic model assumed. 

The systematic error arising from our lack of knowledge of the energy distribution of the events has been estimated 
by repeating the fitting procedure with different energy spectra. From this analysis we estimate a systematic error of 
12%. An additional source of systematic error is related to uncertainties in the hadronic interaction model: following 
the discussion in [1] which relates to the use of QGSJET98 rather than QGSJEtOI, the systematic shift in the fraction 
of iron is 14%. Showers that are calculated using QGSJEtOI are found to develop higher in the atmosphere so that 
the fraction of Fe estimated is reduced from (89 ± 5)% to (75 ± 5)%. 




E(EeV) E(EeV) 

FIG. 8. Fe fraction from various experiments: Fly's Eye (A), Agasa AlOO (■), Agasa Al (□) using SIBYLL 1.5 ([4] and 
references therein) and Haverah Park [1], using QGSJET98 (•)(left) and QGSJEtOI (o)(right). Mean composition determined in 
this paper witli tlie corresponding error for the Volcano Ranch energy range using QGSJET98 (solid line rectangle at left) and 
an estimation of what it would result using QGSJEtOI following [1] (dashed rectangle at right) is shown. 
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VI. COMPARISON WITH OTHER DATA 



In Figure 8 we present the fraction of iron reported by Fly's Eye, Agasa AlOO and Al (using SIBYLL as the hadronic 
model) [6] and the results for Haverah Park [1] using qgsjet98 and qgsjetOI. Also shown is the mean composition 
with the corresponding error we get for the Volcano Ranch energy range using QGSJET98 and an estimation of what 
wc would get using QGSJEtOI as it was done in Haverah Park [1]. 

The inconsistency that exists between several experiments which span different energy ranges, use different tech- 
niques is doubtless enhanced as different hadronic models are used for the interpretation of the raw data. The 
application of a consistent hadronic model [6] brings the results of AGASA into better agreement with the Fly's 
Eye conclusion, while their original analysis stated that there is no indication of changing composition. In Figure 8 
the iron fraction for Fly's Eye and AGASA corresponds to the re-interpretation under the SIBYLL hadronic model, 
including triggering efficiency effects; the original analysis from AGASA is not presented here. 

While Haverah Park, Volcano Ranch and Akeno-AGASA infer X^ax, and hence the overall composition, from 
properties of secondary particle distribution at ground. Fly's Eye and HiRes experiments observe a image of the 
longitudinal shower profile and derive X^ax directly. Nonetheless the estimates can be biased due to the poor 
knowledge of atmospheric properties as recent studies of atmospheric profiles have suggested [28] . 

VII. CONCLUSIONS 

Measurements of the steepness of the lateral distribution r? were made at Volcano Ranch on a shower-to-shower 
basis for fixed bins of zenith angle. We have compared the measured distribution of 77 to our Monte Carlo results for 
proton and iron primaries using QGSJET98 including the scintillator response of the detectors in the Volcano Ranch 
array. Our ability to reproduce Volcano Ranch lateral distribution measurements give us confidence that our analysis 
procedure is correct. 

The cosmic ray mass composition, deduced from Volcano Ranch data, is compatible with mean fraction (89 ± 
5(stat) ± 12(sys)) % of iron in a bi-modal proton and iron mix, in the whole energy range lO^'''^ eV to 10^^ eV, mean 
energy 10^^ eV. Following the discussion in [1], we estimate that this fraction would be reduced to 75 %, with the 
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same qgsjetOI model adopted. 

As mentioned above, some experiments find indications that the composition is becoming lighter with energy. 
However, in a recent work using Haverah Park data above 10^^ eV, it was claimed that the observed time structure 
of the shower front can be best understood if iron primaries are dominant here [27]. The differences between 
measurements on mass composition needs to be addressed further if more solid conclusions on the origin, acceleration 
or propagation of cosmic rays are to be reached. The rate of change with energy and the average mass are still under 
debate. 
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